
#include "gsl_ext_matrix.h"

#include <gsl/gsl_math.h>
#include <math.h>

double gsl_ext_matrix_get_max_abs_diag_element(const gsl_matrix *A)
{
  double max_elem = gsl_matrix_get(A, 0, 0);
  double e;
  int i;
  
  for(i = 1; i < A->size1; i++)
    {
      e = fabs(gsl_matrix_get(A, i, i));
      if(e > max_elem)
	max_elem = e;
    }
  
  return max_elem;
}

double gsl_ext_matrix_get_max_abs_offdiag_element(const gsl_matrix *A)
{
  int i, j;
  double max_elem = gsl_matrix_get(A, 0, 1);
  double e;
  
  for(i = 0; i < A->size1; i++)
    {
      for(j = 0; j < A->size2; j++)
	{
	  if(i != j)
	    {
	      e = fabs(gsl_matrix_get(A, i, i));
	      if(e > max_elem)
		max_elem = e;
	    }
	}
    }
  
  return max_elem;
}

void gsl_ext_matrix_compute_gerschgorin_bounds(const gsl_matrix *A, double *minev, double *maxev)
{
  int i, j;
  double row_sum;
  *minev = *maxev = gsl_matrix_get(A, 0, 0);
  double e;
  
  for(i = 0; i < A->size1; i++)
    {
      row_sum = 0.0;
      for(j = 0; j < i; j++)
	row_sum += fabs(gsl_matrix_get(A, i, j));
      for(j = i + 1; j < A->size2; j++)
	row_sum += fabs(gsl_matrix_get(A, i, j));
    
      e = gsl_matrix_get(A, i, i);
      *maxev = GSL_MAX_DBL(*maxev, e + row_sum);
      *minev = GSL_MIN_DBL(*minev, e - row_sum);
    }
}
